## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Attachement du package : 'plotly'
## 
## 
## L'objet suivant est masqué depuis 'package:ggplot2':
## 
##     last_plot
## 
## 
## L'objet suivant est masqué depuis 'package:stats':
## 
##     filter
## 
## 
## L'objet suivant est masqué depuis 'package:graphics':
## 
##     layout
main_theme = theme_bw()+
  theme(line = element_blank(), 
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.ticks =  element_line(colour = "black"),
        axis.text.x = element_text(colour = "black", size=22),
        axis.text.y = element_text(colour = "black", size=22),
        legend.title = element_text(colour = "black", size=20),
        legend.title.align=0.5,
        legend.text = element_text(colour = "black", size=18),
        axis.title=element_text(size=28))
#Parametres
i = 20 # Nombre d'especes
Xmin = 0
Xmax = 2
X0 = (Xmax - Xmin)/2
K0 = 1
lambda = 1
sigma = 1 #etalement de la competition
r = 1
############ Fonctions #############

K_cap = function(x, K0_ = K0, lambda_= lambda, x0 = X0){
   (K0_ - lambda_*(x-x0)**2)
}

a = function(x1,x2, sigma_a = sigma){ 
  exp( (-1/2*(x2-x1)**2)/sigma_a**2 )
}

fitness = function(x1, x2, r_ = r ){
  r_*(1-a(x1, x2)*(K_cap(x1)/K_cap(x2)))
}
## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 25);
x2_interval  <- seq(0.1, 1.9, length.out = 25);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

p  <- persp(x1_interval,x2_interval,z, theta = 30, phi = 20,
            col = "lightblue", shade = 0.4, ticktype = "detailed")

p
##               [,1]       [,2]       [,3]       [,4]
## [1,]  9.622504e-01 -0.1900112  0.5220515 -0.5220515
## [2,]  5.555556e-01  0.3291090 -0.9042196  0.9042196
## [3,] -1.567806e-17  0.4812177  0.1751489 -0.1751489
## [4,] -1.517806e+00  0.3651056 -2.1663676  3.1663676
plot_ly() %>% add_surface(x = ~x1_interval, y = ~x2_interval, z = ~z)
## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 500);
x2_interval  <- seq(0.1, 1.9, length.out = 500);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

couleurs <- ifelse(z > 0, "lightblue", "black")

matrice_couleurs = couleurs
df <- data.frame(x = rep(1:nrow(matrice_couleurs), each = ncol(matrice_couleurs)),
                 y = rep(1:ncol(matrice_couleurs), times = nrow(matrice_couleurs)),
                 couleur = as.vector(matrice_couleurs))

ggplot(df, aes(x = x, y = y, fill = couleur)) +
  geom_tile() +
  scale_fill_identity() +
  labs(title = "PIP") +
  theme_minimal()+
  labs(fill = "Légende\nNoir = Négatif\nBleu = Positif")

K_cap = function (X){
 max(0, K0 - lambda*(X - X0)^2) + 10^(-9) # On rajoute 10^-9 pour ne pas avoir de pb de division par 0
}
sigma = 0.3
############### Simuler LV pour un nombre arbitraire d'espece ############

#Conditions initiales
x = seq(from = 0, to = 2, length.out = i) # Valeurs des phenotypes
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
 # K de chaque espece
n = rep(1/i,i) #densite de pop initiale

# On met sous forme de matrice
N0 = t(as.matrix(n))
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
X = matrix(rep(x,i), ncol = i, nrow = i, byrow = TRUE)

#Calcul de la matrice M
D = t(X) - X
A = exp(-0.5*D^2/(sigma^2))
M = A/t(K)

#Resolution du systeme d'ODE
model <- function(t,N,sigma){
  dN <- r * N * (1 - N %*% M)
  return(list(dN))
}

ode <-ode(N0,0:100,model,parms = sigma)
# Representation
#barplot(ode[100,1:i+1])

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))

ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]

ggplot(ode_plot)+
  geom_raster(aes(trait, time,  fill=dens))+
  scale_fill_gradient2(low = "white" ,
                       high = "red")

######## Nombre d'espece qui persiste en fonction de sigma (pour l'instant ne marche pas..) ##############
seuil <- 0.05
nb_especes <- c()
for (sigma in seq(0.3,5,by = 0.1)){
  ode <-ode(N0,0:100,model,parms = sigma)
  nb = 0
  for (m in tail(ode,1)[,1:i+1]){
    if (m>seuil){
      nb = nb + 1
    }
  }
  nb_especes <- c(c(nb_especes),c(nb))
  
}
plot(seq(0.3,5,by = 0.1),nb_especes)

######## Implementation de mutation ##########

#parametres
p = 0.1 #taux de la pop qui mute
T = 10000
t = 200
n = rep(0,i)
n[i/6] <- 1 # densites de pop initiale (une seule espece)
N0 <- t(as.matrix(n))

#Initialisation des valeures de K
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)

M = A/t(K)

ode <-ode(N0,1:t,model,parms = sigma)
ode_mutation <- ode

#trouver un moyen de faire muter tout le monde... (peut etre metre un seuil de densité min..)

for (z in 2:(T/t)){
  N0 <- t(as.matrix(tail(ode,1)[,1:i+1]))
  a = runif(1)
  index = which.max(N0)
  if (a<0.5){ #mutation vers la gauche
    N0[max(0,index-1)] <- p*N0[index]
    N0[index] <- (1-p)*N0[index]
  } else { #mutation vers la droite
    N0[min(i,index+1)] <- p*N0[index]
    N0[index] <- (1-p)*N0[index]
  }
  ode <-ode(N0,(z*t):(z*t+t),model,parms = sigma)
  
  ode_mutation <- rbind(ode_mutation,ode)
}

ode_mutation_plot = as.data.frame(ode_mutation)%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ode_mutation_plot$sp_ID <- as.numeric(ode_mutation_plot$sp_ID)
ode_mutation_plot$trait <- x[as.numeric(ode_mutation_plot$sp_ID)]

ggplot(ode_mutation_plot)+
  geom_raster(aes(trait, time,  fill=dens))+
  scale_fill_gradient2(low = "white" ,
                       high = "red") +
  main_theme
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

# ggplot(ode_mutation_plot)+
#   geom_tile(aes(trait, time,  fill=dens))+
#   scale_fill_gradient2(low = "white" ,
#                        high = "red") +
#   main_theme